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Abstract. 

We simulate the dynamics of fractal star clusters, in order to investigate the evolution of substructure in recently formed clusters. 
The velocity dispersion is found to be the key parameter determining the survival of substructure. In clusters with a low initial 
velocity dispersion, the ensuing collapse of the cluster tends to erase substructure, although some substructure may persist 
beyond the collapse phase. In clusters with virial ratios of 0.5 or higher, initial density substructure survives for several crossing 
times, in virtually all cases. Even an initially homogeneous cluster can develop substructure, if it is born with coherent velocity 
dispersion. 

These results suggest that the simple initial conditions used for many sophisticated A'-body simulations could be missing a very 
important and dramatic phase of star cluster evolution. 
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1. Introduction 

It appears that most stars - possibly all stars - form in clus- 
ters. Their dynamical evolution is thus of great interest. In re- 
cent years, codes such as nbody6 (Aarseth 2000), which include 
detailed stellar evolution and mass loss, binary evolution and 
mass transfer, have made possible a new generation of 'kitchen 
sink' simulations (e.g. Kroupa et al. 2001; Hurley et al. 2001; 
Portegies Zwart et al. 2001). However, the initial conditions 
for these simulations have often been very simple, for example 
Plummer models, in stark contrast with the great detail invoked 
in modelling the subsequent evolution. 

Observations of star clusters, on the other hand, suggest that 
the initial conditions for star formation are highly clumpy and 
structured, both in the distribution of the molecular gas from 
which stars are about to form (e.g. Williams 1999 and refer- 
ences therein), and in the distribution of newly-formed stars 
(e.g. Bate et al. 1998; Gladwin et al. 1999). 

Aarseth & Hills (1972) were the first to investigate the evo- 
lution of collapsing star clusters with substructure. Their sim- 
ulations were limited by the available computer power to clus- 
ters of 120 stars. They found that subclustering was destroyed 
on a free-fall timescale. Later Goodwin (1998) investigated an 
initially virialised cluster with density substructure and a larger 
number of stars. He found that most of the initial substructure 
was erased within a few crossing times. 

In this paper we investigate the evolution of initially frac- 
tal star clusters to see how long substructure can survive. Our 
models include star clusters with large velocity dispersions, 
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as would be expected in clusters shortly after they expel their 
residual gas (cf. Goodwin 1997). In addition, we investigate 
fractal clusters in which the density substructure is correlated 
with coherent velocity dispersion, as would be expected in clus- 
ters where sub-clusters form from distinct molecular cores. 

We are not suggesting that star clusters are necessarily frac- 
tal. The range of scales over which young star clusters exhibit 
substructure is usually very small, often less than an order of 
magnitude, so the notion of a fractal cannot be applied rigor- 
ously. Nonetheless, fractals provide a simple, one-parameter 
description of dumpiness, and this is why we are using them. 

2. N-body method and initial conditions 

We conduct our A^-body simulations on a grape-5 board, which 
allows for very rapid solution of the A^-body gravitational prob- 
lem (Kawai et al. 2000). We use a simple direct first order N- 
body integrator, as the speed of the grape-5 board allows the 
timestep to be set sufficiently small that over the course of a 
simulations the total energy of the system never changes by 
more than 0.01 per cent, and usually by significantly less than 
this. A small softening length is used (generally of order 10~ 5 
in code units). 

The presence of binaries in a cluster will effect the dy- 
namical evolution of the system by altering impact parameters. 
However, tests using different softening lengths show the re- 
sults to be independent of the softenning length, even when it 
is significantly larger than our canonical 10~ 5 . This is due to 
the relaxation being primerally violent relaxation, rather than 
encounter-driven two-body relaxation. 
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Fractal distributions are generated by defining an Mr-cube 
with side 2, and placing an Mr-parent at the centre of the ur- 
cube. Next, the ur-cube is divided into iVl equal sub-cubes, 
and a child is placed at the centre of each sub-cube (the first 
generation). Normally we use N,a v = 2, in which case there 
are 8 sub-cubes and 8 first-generation children. The probabil- 
ity that a child matures to become a parent in its own right 
is M D 3) , where D is the fractal dimension; for lower D, the 

div ' ' ' 

probability that a child matures to become a parent is lower. 
Children that do not mature are deleted, along with the ur- 
parent. A little noise is then added to the positions of the 
remaining children, to avoid an obviously gridded structure, 
and they become the parents of the next generation, each one 
spawning children (the second generation) at the centres of 
A^j equal-volume sub-sub-cubes, and each second-generation 
child having a probability N ( £~ 3) of maturing to become a par- 
ent. This process is repeated recursively until there is a suf- 
ficiently large generation that, even after pruning to impose 
a spherical envelope of radius 1 within the ur-cube, there are 
more children than the required number of stars. Children are 
then culled randomly until the required number is left, and the 
survivng children are identified with the stars of the cluster. 

We explore a range of models with N tot = 1000 and 
N to t = 10000 stars. The fractal dimensions which we investi- 
gate are D = 1.6, 2.0, 2.6, &3.0, since these all correspond 
to 2 D (the mean number of maturing children) being close to 
an integer. This reduces the likelihood of departures from the 
specified fractal dimension, because in our algorithm for con- 
structing an initial cluster the number of maturing children born 
to each parent should be an integer. 

Our simulations begin with a random velocity dispersion, 
which is either incoherent or coherent. For an incoherent veloc- 
ity dispersion, each particle is given random cartesian velocity 
components from a gaussian distribution, and these velocities 
are then scaled so that the virial ratio a has the prescribed value. 
a is defined as the ratio of the total kinetic energy to the magni- 
tude of the gravitational potential energy, a = 0.1 corresponds 
to a cluster which immediately collapses, a = 0.5 corresponds 
to virial equilibrium, a = 0.75 corresponds to a super-virial 
cluster (i.e. one which has just expelled its residual gas). 

For a coherent velocity dispersion, each star inherits most 
of its velocity from its family tree. We first calculate, for each 
mature child of each generation, how many stars are descended 
from it, and this gives the child's mass. Next, starting with the 
children of the first generation, we give them random velocities 
relative to the ur-parent, so that they form a virialised cluster. 
We treat the children of each subsequent family in the same 
way, giving them random velocities relative to the parent, so 
that they form a virialised sub-cluster, but in addition they ac- 
quire the velocity of the parent. This process is repeated recur- 
sively, and hence the stars of the final generation inherit random 
velocities from all their antecedents, giving a coherent velocity 
dispersion. Finally, the velocities are scaled so that the virial 
ratio a has the prescribed value. 

We believe that the most realistic of the initial velocity dis- 
persions we use is the coherent, super-virial one (a = 0.75). 
Our simulations do not include a gaseous component, but we 



assume that the residual gas has been expelled from the cluster 
very rapidly and recently. If the stars have previously been in 
virial equilibrium with the residual gas, the virial ratio of the 
purely stellar component is super-virial. It also seems likely 
that the velocities of the stars will be correlated locally, on the 
assumption that each sub-cluster of stars has formed from a 
single molecular core. 

3. Measures of dumpiness 

The primary aim of this paper is to investigate how rapidly 
initial density substructure and associated velocity coherence 
are destroyed, and the typical timescale before a cluster can be 
described as smooth. In order to do this we require a robust 
method for measuring the dumpiness of a cluster. 

One method would be to determine the evolution of the 
fractal dimension of the cluster. However, this method has 
drawbacks, because it is difficult to measure the fractal dimen- 
sion of a distribution with a strong overall radial density gradi- 
ent. Clusters rapidly acquire a core-halo structure, in which the 
mean separation between stars decreases radially by a signifi- 
cant factor. To determine the fractal box-dimension of a cluster, 
a grid of cells is overlaid on the cluster, and the number of cells 
containing at least one star is counted. This is repeated, start- 
ing with a coarse grid of large cells, and proceeding to ever 
finer grids of smaller cells. The fractal box-dimension is then 
given by the slope of a plot of the log of the number of oc- 
cupied cells against the log of the inverse cell size. This slope 
turns over at the point where the grid becomes saturated, i.e. 
where the number of occupied cells is equal to the number of 
stars; further reductions of cell size then make no difference to 
cell occupancy. However, in the presence of a density gradient, 
saturation of the halo occurs before saturation of the core, pro- 
ducing a complicated relationship that is difficult to interpret. 
It might be possible to compensate for this, if the overall radial 
density distribution is known a priori, but in general it is not. 

The simple method we use to measure dumpiness is to ac- 
cord each star a local density, given by Sm/Vs, where m is the 
stellar mass and Vs is the spherical volume bounded by the 
fifth nearest star. Next we overlay spherical shells centred on 
the centre of mass of the ten densest stars, and work out the av- 
erage density in each shell. Substructure is then visible as stars 
(or groups of stars) within a shell having density significantly 
higher than the average density for that shell. As a measure of 
the level of dumpiness, we use two numbers, F20 and F50. F20 
(F50) measures the fraction of shells in which more than 20% 
(50%) of the stars have more than 5 times the average density. 
A high level of dumpiness is reflected in large values of F20 
and F 5 q. We tested more complicated, kernel weighted density 
estimates (such as those used in SPH simulations), but found 
no significant difference in the results. 

F20 and F50 are not greatly affected by small number statis- 
tics. Tests performed with clusters of 1000 stars show that 
F50 < ^20 ~ 0.1 for various clusters which should not have 
significant dumpiness, viz. a randomly distributed uniform- 
density cluster, a randomly distributed Plummer cluster, and 
a D = 3 fractal cluster. In contrast, a D = 1.6 fractal clus- 
ter, which should be very clumpy, has F50 ~ F20 ~ 1- We 
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conclude that F20 and F50 are useful measures of dumpiness. 
Figure[Jshows a fractal distribution with obvious substructure, 
and the large points show which stars have been selected by 
this method as being 5 times the average density. They coin- 
cide well with the substructure which the human eye identifies. 

4. The evolution of fractal clusters 

We now use F20 and F50 to investigate the evolution and erasure 
of substructure from an initially fractal cluster. 

4.1. Collapsing clusters 

Figure |2] shows typical results for the evolution of F20 and 
F50 for clusters having initial virial ratio a — 0.1 (such 
that they will collapse) and initial fractal dimensions D = 
1.6, 2.0, 2.6, & 3 (decreasing levels of initial substructure). In 
all of these cases the initial velocity dispersion is coherent. 

Time is given in A^-body units, such that G — M = R = 
1, where M and R are the total mass and initial radius of the 
cluster. For example, if the total mass of the cluster is M = 
100M o and the initial radius is R — 1 pc, then one time unit is 
« 1.5 Myr (see Heggie & Mathieu 1986 for details of N-body 
units). 

Figure |2] shows that the level of substructure tends to de- 
crease with time, and all substructure has essentially disap- 
peared by T = 3. An interesting feature is the transient rise in 
the level of substructure in clusters with high D (i.e. those with 
initially low levels of substructure). This is due to the coherence 
in the initial velocity field. Even though there is initially a low 
level of density substructure, for a short time the coherent ve- 
locity field generates substructure. However, this substructure 
does not last long. 

The main mechanism for erasing substructure is the gravi- 
tational interactions between clumps. The potential of a clumpy 
cluster is highly uneven and violent relaxation occurs, allowing 
the stars to relax into a smooth distribution on a short timescale. 
Two-body encounters also act to remove kinetic energy from 
the main cluster by ejecting stars. This can be seen in a rapidly 
expanding halo of stars around the main cluster, and a cluster 
core which is often significantly smaller than the initial size of 
the system. 

Figure [5] shows in detail the evolution of the D — 2 clus- 
ter from the previous figure. By inspection the evolution of the 
substructure is well described by our measures F20 and F50. In 
Fig.[3]the initially very clumpy distribution rapidly collapses. 
At first some of the clumps disperse, but those which are bound 
collapse and the larger of these clumps then attract nearby 
clumps. The initial dumpiness is erased as most clumps merge. 
However, one of the merging clumps is sufficiently bound to 
survive the initial merger process and emerges at late times to 
the bottom left. The re-emergence of this clump explains the 
increase in F20 and F50 for this cluster in Fig.|5]around T = 2. 

4.2. Virialised clusters 

Figure @] shows the evolution of F20 and F50, as in Fig. [5] but 
for clusters with an initial virial ratio of a = 0.5. Clusters with 



a = 0.5 have enough kinetic energy to support themselves 
against overall collapse, but their clumpy nature means that sig- 
nificant dynamical evolution still occurs. The most significant 
features of Fig.0]are the longer time required for substructure 
to be erased, and the persistence of substructure throughout the 
simulation in the D — 1 .6 cluster. 

The initial coherence of the velocity dispersion in the high- 
D simulations is again responsible for the increase in sub- 
structure early on, although this substructure is erased over a 
timescale of a few time units (a few Myr in a typical cluster). 

In the D = 1 .6 cluster the survival of substructure is due to 
the high velocity dispersion and its coherence. As can be seen 
in Fig. [5] the initial cluster rapidly divides into a main cluster 
and a lumpy sub-cluster, which has sufficient bulk velocity to 
escape from the main sub-cluster, dividing into 3 small sub- 
sub-clusters as it does so. 

Such an evolution is typical for clusters of low fractal di- 
mension and not unusual for clusters of higher fractal dimen- 
sion. The key is the coherent nature of the velocity dispersion. 
The velocity dispersion within a sub-cluster holds it up against 
collapse, and its bulk velocity helps it to avoid merging with 
other sub-clusters. Hence it is able to maintain its separate iden- 
tity for a long time. This contrasts with the cases discussed in 
the previous subsection, where the collapse tends to erase sub- 
structure rather quickly. 

4.3. Supervirial clusters 

Figure. [6] again shows F20 and F50, this time for clusters with 
an initial virial ratio of a — 0.75, such that they expand. In 
these cases it is clear that the level of substructure does not 
decrease rapidly. Significant substructure remains at the end of 
the simulations for all D. 

Even clusters whose density structure is initially not very 
clumpy grow density substructure, again due to the coherence 
of their initial velocity dispersion. Generally, the lower D, the 
more clumps there are, and the smaller they are. For high D 
often a binary cluster is formed with two major sub-clusters. 
Figure shows the final states (at T = 5) for different realisa- 
tions of clusters with a variety of initial fractal dimensions. In 
virtually all simulations clusters still have a significant amount 
of well defined substructure present at T = 5. (For total mass 
M « IOOOMq and initial radius R a 1 pc, T — 5 corresponds to 
5 to 10 Myr.) 

Typically, by T — 5 a main sub-cluster is identifiable (as 
a sub-cluster that is significantly larger than all other sub- 
clusters). Figure0shows in most cases a clear main sub-cluster 
with surrounding substructure. In many cases the substructure 
is never erased altogether; some clumps are not bound to the 
main clump and so escape. This is very common in clusters 
with an initial virial ratio of a — 0.75, and even occurs oc- 
casionally in clusters with an initial virial ratio of a = 0.5. 
Sub-clusters that are bound to the main sub-cluster can remain 
as separate entities for a significant length of time. We ran a 
subset of our simulations until T — 100 and found some sub- 
clusters remaining in orbit. However, in the majority of cases 
the tidal effect of the main sub-cluster disrupts any bound sub- 
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structure by T — 10 or 20 (roughly 10 to 50 Myr). This leaves 
a smooth, but typically elliptical, main cluster surrounded by 
a large, and expanding, halo of stars. We did not however in- 
vestigate the details of this phase of cluster evolution in any 
detail. 

4.4. Incoherent velocity structures 

Given the potential importance of coherent velocity dispersion 
in maintaining - and even increasing - the level of substructure, 
it is expected that clusters with incoherent velocity dispersion 
will not retain substructure so long. Figure [8] shows the evolu- 
tion of F20 and f 50 for a cluster with initial virial ratio a = 0.75 
and incoherent velocity dispersion. The substructure is erased 
almost immediately. This is unsurprising, as, with such a large 
virial ratio, any density substructure almost instantly disperses. 

However, we believe that coherent velocity dispersion is the 
more plausible initial condition. The likely cause of the sub- 
structure that is commonly observed in young star clusters is 
the formation of sub-clusters of stars in distinct molecular cores 
formed by the turbulence in the molecular cloud. It is therefore 
to be expected that the stars within a sub-cluster have similar 
velocities, and that the velocity dispersion is correlated with the 
density substructure. 

4.5. The effect of large N 

The previous results have concentrated on clusters with N tot = 
1000 representing moderately rich, open clusters (for example, 
the Pleiades or Orion). Here we investigate the effect of larger 
N tot on the survival or destruction of substructure. 

Figure [9] shows the evolution of F20 and F50 for clusters 
with virial ratio a = 0.75, coherent velocity dispersion, and 
D = 1.6, 2, 2.6, &3, equivalent to Fig. ©but with AT tot = 
10000. Figure [5] illustrates a common theme for large N tot , 
namely that the decline in F20 and F50 is somewhat more rapid 
than in clusters with small N tot . 

This result is largely an artifact of the clump detection pro- 
ceedure. Often in N - 1000 clusters the number of particles 
within a particular radial bin is small (of order a few), whilst in 
N = 10000 clusters that number is (obviously) far larger. Thus 
poisson noise has more of an effect upon the statistics in low-N 
clusters. 

As an example, two a = 0.75, D — 2 clusters were simu- 
lated, one with N = 10000 and the other with 90 per cent of 
the particles from the first simulation removed at random. At 
T = 5, the N = 1000 cluster had F 20 = 0.85 and F 50 = 0.25, 
and for N = 10000, F 20 = 0.70 and F 50 = 0.05. The F- 
statistics in the N = 1000 simulation were increased by the 
effect of low numbers of particles (5 to 10) in the inner radial 
bins, while the N = 10000 cluster suffered far less from this 
effect with no less than 20 particles in any one bin. The den- 
sity of particles in low-N simulations is also less smooth than 
in high-Af simulations as the search radius in which 5 particles 
are to be found is usually significantly greater. Small clumps in 
low-N simulations tend to be more above the average density 
than in large-Af simulations. 



Nonetheless, large N tot clusters still maintain significant 
levels of substructure for several crossing times. This may be an 
explanation for some of the anomalous bumps observed in the 
profiles of several LMC star clusters (e.g. Mackey & Gilmore 
2003). 

4.6. The final structure of clusters 

An examination of the most significant sub-clusters at the end 
of simulations shows that they generally appear similar to older 
clusters. The relaxation of the clusters creates Plummer- or 
King-like profiles, ie. a flat central density profile, with a steep 
decline in the halo. This is especially clear in clusters with 
N = 10000. The cluster illustrated in Fig. ??, has a flat cen- 
tral surface density of ss 500 stars per unit area over a radius of 
~ 0.5, followed by an approximate r~ 2 decline. There are sig- 
nificant sub-clusters apparent in the halo, the largest of which 
cause 'bumps' to appear in the surface density profile. For clus- 
ters with low N, the trend to form core-halo density structures 
is present, but far less clear due to low -AT noise. 

Some clusters have distinct, unbound sub-clusters which 
are included in our dumpiness determinations. Once a sub- 
cluster has travelled a significant distance from the main cluster 
it will be difficult (without proper motions) to determine that it 
formed in the same position as the main cluster. Thus it may 
appear as if they are smooth and separate clusters that formed 
coevally, rather than the result of dynamical segregation from 
the same initial cluster. 

5. Implications for observations and simulations 

Observations of young star clusters often show a very inho- 
mogeneous, clumpy distribution. Our simulations demonstrate 
that if the velocity dispersion of these clusters is low, then much 
of that initial substructure will be erased in the ensuing collapse 
(cf. Aarseth & Hills 1972). However, if, as we might expect, the 
velocity dispersion is high, such that the cluster remains sup- 
ported or even expands, then significant levels of substructure 
can survive for several crossing times. Even an initially homo- 
geneous distribution of stars can grow substructure if the initial 
velocity dispersion is coherent. 

These results have two important consequences. First, 
young clusters probably undergo significant and rapid dynam- 
ical evolution. Therefore drawing conclusions about the ini- 
tial stellar distribution from observations is very difficult, and 
should take these considerations into account. Taurus is a good 
example of a young, embedded cluster with a high level of sub- 
structure (Briceno et al. 1993; Ghez 1993). Our results imply 
that, when Taurus expels its residual gas and becomes a pure N- 
body system, it will probably not simply collapse into a small, 
dense open cluster. Indeed, it is more likely that it will separate 
into two or three small clusters, which in a few Myrs may look 
as though they formed separately from each other. As we have 
shown, this depends crucially on the virial ratio of the final (gas 
free) cluster, and on the coherence of the velocity dispersion. 
Even a cluster which is initally fairly smooth in appearance, 
such as IC348 (Najita et al. 2000), may develop substructure 
and look very different in a few Myrs. 
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Second, the rapid dynamical evolution early in a cluster's 
life is also important when setting the initial conditions for N- 
body simulations of cluster evolution. Simple initial conditions 
(such as Plummer spheres) will fail to capture this potentially 
important stage of cluster evolution. Processes like binary cap- 
ture or dissolution, and stellar ejection, which probably occur 
in the first few Myr of a cluster's lifetime, may be strongly af- 
fected by the evolving density substructure, and the dynamical 
changes we have identified here. 
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Fig. 1. Left: a clumpy distribution of stars. Right: the stars selected by our dumpiness measure as being significantly overdense. 
This distribution has high F20 and Fso which indicate that it is significantly clumpy (which the eye confirms). 
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Fig. 2. The evolution of the F20 and F$q measures of dumpiness for clusters with 1000 particles, initial virial ratio 0.1, coherent 
velocities and initial D of 1.6 (solid lines), 2 (dashed lines), 2.6 (dot-dashed lines) and 3 (dotted lines). 



Simon P. Goodwin and A. P. Whitworth: The dynamical evolution of fractal star clusters: the survival of substructure 7 









Fig. 3. The evolution of an initially Fdim = 2 cluster with N = 1000 and a virial ratio a = 0.1 with coherent velocities. The time 
in iV-body units is given in the top right of each panel. 
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Fig. 5. The same as Fig.|3]but for a cluster with an initial virial ratio of a = 0.5. 
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Fig. 7. The final states of different realisations of clusters at T = 5, with N = 1000, virial ratios of a = 0.75 and coherent 
velocities. The initial fractal dimension is labeled in each case in the top right. 
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8. As Fig.|5]but for clusters with an initial virial ratio of a — 0.75 and a random rather than coherent velocity substructure. 
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9. As Fig.|6]but for clusters with N = 10000. 



